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IDENTIFICATION  OF  CONNECTIVITY 
IN  NEURAL  NETWORKS 


Xiaowei  Yang  and  Shihab  A.  Shamma 
Electrical  Engineering  Department,  Systems  Research  Center, 
and  the  University  of  Maryland  Institute  for  Advanced  Computer  Studies 
University  of  Maryland,  College  Park,  Maryland  20742,  USA 


Abstract.  Analytical  and  experimental  methods  are  provided  for  estimating 
synaptic  connectivities  from  simultaneous  recordings  of  multiple  neurons.  The  results 
are  based  on  detailed,  yet  flexible  neuron  models  in  which  spike  trains  are  modeled  as 
general  doubly  stochastic  point  processes.  The  expressions  derived  can  be  used  with 
non-stationary  or  stationary  records,  and  can  be  readily  extended  from  pair-wise  to 
multi-neuron  estimates.  Furthermore,  we  show  analytically  how  the  estimates  are 
improved  as  more  neurons  are  sampled,  and  derive  the  appropriate  normalizations  to 
eliminate  stimulus-related  correlations.  Finally,  we  illustrate  the  use  and  interpreta¬ 
tion  of  the  analytical  expressions  on  simulated  spike  trains  and  neural  networks,  and 
give  explicit  confidence  measures  on  the  estimates. 


1 


1  Introduction 


Most  functions  of  the  mammalian  nervous  system  are  performed  by  networks 
of  highly  interconnected  neurons.  In  the  experimental  study  of  these  networks,  ex¬ 
tracellular  recording  are  often  employed  to  sample  the  patterns  of  action  potentials 
simultaneously  generated  by  several  neurons  (Abeles  and  Goldstein  1977;  Gerstein 
and  Perkel  1969;  Yang  and  Shamma  1988).  The  correlations  among  the  recorded  fir¬ 
ings  of  the  different  cells  are  then  used  as  measures  of  the  type  and  strength  of  their 
interconnections.  Many  such  measures  have  been  proposed  to  accomplish  the  latter 
task;  they  include  the  cross-interval  histograms,  the  cross-correlation  histograms,  the 
cross-covariance  histogram,  and  the  joint  peri  stimulus  time  (PST)  histogram  (the 
scatter  diagram)  (Gerstein  1970;  Gerstein  and  Perkel  1969).  In  all  cases,  the  his¬ 
tograms  provide  statistical  measures  in  support  of  various  hypotheses  such  as  whether 
the  two  (or  more)  neurons  under  study  directly  influence  each  other  or  simply  share 
common  inputs,  and  whether  the  influences  are  excitatory  or  inhibitory. 

There  are  three  basic  difficulties  with  these  methods  that  we  tackle  in  this  re¬ 
port.  The  first  concerns  the  lack  of  flexible  general  analytical  treatments  that  outline 
the  relations  between  the  synaptic  connectivities  and  the  correlation  measures  that 
are  used  to  estimate  them.  Thus,  while  various  features  in  the  above  mentioned 
histograms  may  reflect  qualitatively  the  underlying  connections,  several  parameters 
and  conditions  can  render  these  measures  inadequate.  Examples  of  such  difficulties 
are  the  differing  integrating  dynamics  of  different  cell  types,  and  the  potentially  se¬ 
vere  errors  due  to  stimulus-induced  (rather  than  synaptic)  correlations.  Attempts  to 
overcome  these  problems,  as  in  the  use  of  the  shuffling  method  to  reduce  stimulus 
effects,  are  shown  here  to  be  largely  inadequate. 

The  second  basic  shortcoming  of  the  above  correlation  methods  stems  from 
the  nonstationarity  of  the  neural  records.  In  constructing  cross-interval  and  cross- 
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correlation  histograms,  counts  are  usually  obtained  not  only  by  averaging  over  dif¬ 
ferent  stimulus  presentation  but  also  by  averaging  over  the  time  duration  of  each 
presentation  period.  This  makes  these  two  estimates  inadequate  when  working  with 
non-stationary  records  and,  instead,  measures  based  on  time-dependent  histograms 
such  as  the  joint  PST  scatter  diagram  should  be  used  for  the  analysis  (Habib  and 
Sen  1985;  Stokkum  et  al.  1986). 

Finally,  it  is  unclear  in  many  existing  methods  how  to  extend  the  analysis  to  more 
than  two  neurons,  and  how  to  evaluate  the  degree  to  which  a  pair-wise  estimate  is 
improved  when  the  records  from  many  other  neurons  are  included.  This  is  a  partic¬ 
ularly  important  criterion  as  progress  in  multiunit  recording  technologies  promises 
to  increase  significantly  the  number  of  records  of  simultaneously  active  neurons. 

To  summarize,  the  objectives  of  this  report  are  (1)  to  provide  rigorous  analyti¬ 
cal  and  experimental  methods  to  estimate  synaptic  connectivities  from  simultaneous 
recordings  of  multiple  neurons  that  are  based  on  accurate  and  flexible  neuron  mod¬ 
els;  (2)  to  express  synaptic  connectivity  in  terms  of  probability  densities  of  joint 
neuronal  firings  and  individual  neuronal  firings  that  can  be  used  with  non-stationary 
(or  stationary)  records;  (3)  to  extend  these  methods  from  pair-wise  to  multiunit 
correlations. 

The  paper  is  organized  as  follows.  In  the  next  section,  a  stochastic  nonlinear 
neuron  model  is  proposed,  and  the  spike  train  generated  by  the  model  is  expressed 
by  a  doublely  stochastic  process.  This  model  will  serve  as  the  fundamental  tool 
upon  which  the  analytical  results  are  based.  In  section  3,  quantitative  analyses  of 
neuronal  connectivities  are  carried  out  through  the  model.  These  include  derivations 
of  the  relations  between  the  synaptic  connectivity  and  the  firing  probability  densities, 
and  extending  the  pair-wise  correlations  to  the  multi-neuron  case.  In  section  4,  the 
results  are  summarized  and  discussed  in  the  context  of  practical  implementations 
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and  considerations  of  the  accuracy  of  the  estimates.  Finally,  the  analytical  results 
are  simulated  and  discussed  in  section  5.  The  proofs  of  lemmas  and  theorems  are 
given  in  the  Appendix. 

All  the  analytical  treatments  are  contained  within  sections  2  and  3.  For  the 
reader  interested  only  in  using  the  final  expressions,  section  4  outlines  the  results 
and  is  sufficient  as  a  guide  for  their  experimental  applications. 


2  The  Neuron  Model 

The  basic  unit  of  the  nervous  system  which  receives  and  transmits  neural  signals 
is  the  neuron.  The  interactions  of  neurons  in  a  network  occur  in  most  cases  through 
synaptic  connections  between  them.  Most  synapses  are  found  between  the  axon 
terminals  of  a  presynaptic  neuron  and  the  soma  or  dentritic  tree  of  a  post-synaptic 
neuron.  Since  there  can  be  many  synapses  between  any  two  neurons,  it  is  impractical 
in  modeling  the  neural  network  to  account  for  individual  synapses;  rather,  it  is  more 
fruitful  both  for  experimental  investigation  and  mathematical  description  to  consider 
the  total  effective  influence  of  one  cell  on  another. 

A  synaptic  connection  from  a  presynaptic  neuron  B  to  post-synaptic  neuron  A  is 
said  to  be  excitatory  (inhibitory)  if  the  firing  rate  of  neuron  A  increases  (decreases) 
when  neuron  B  fires.  For  the  purposes  of  the  model,  we  assume  that  the  post- 
synaptic  potentials  due  to  many  presynaptic  inputs  are  continuously  integrated  to 
produce  a  post-synaptic  membrane  potential.  A  neuron  fires  an  action  potential  when 
its  membrane  potential  exceeds  a  threshold  level.  After  each  action  potential,  there 
is  a  period  during  which  the  probability  of  firing  is  reduced.  This  period  is  divided 
into  two  intervals:  The  first  is  the  absolute  refractory  period ,  in  which  the  neuron 
cannot  fire  again;  the  second  is  the  relative  refractory  period  where  the  neuron  may 
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generate  a  spike  only  when  the  stimulus  is  fairly  strong. 


Since  the  action  potentials  of  a  given  neuron  are  similar  in  shape,  we  assume  that 
the  transmitted  information  is  carried  only  through  the  temporal  patterns  of  the 
spike  trains,  and  hence  we  use  a  sequence  of  impulses  to  abstract  a  train  of  action 
potentials.  Because  the  instantaneous  firings  of  a  neuron  are  not  deterministic,  a 
stochastic  point  process  is  adopted  to  model  the  firings  (Perkel  et  al.  1967).  All 
the  stochastic  processes  and  random  variables  to  be  discussed  are  defined  on  some 
probability  space  (17,  F,  P).  Let  (0,  F,  P)  be  a  probability  space,  and  let  {Ft  :t>  0} 
be  a  non-decreasing  family  of  sub  cr-field  of  F  (i.e.,  F„  C  Ft,  for  every  s  <  t).  The 
family  {Ft}  is  called  a  history,  and  Ft  represents  the  information  collected  during 
[0,7],  Let  {Vj}t>o  be  a  stochastic  process  (representing  the  semi-membrane  potential 
process)  defined  on  (£l,F,P).  The  family  of  the  sub  cr-field  generated  by  {Vt}t>o, 
7 {t  =  o(Vs  :  0  <  s  <  t),  is  called  a  history  of  {V)}i>o  if  Wt  C  Ft,  for  t  >  0.  And  in 
this  case,  {Vi}t>o  is  said  to  be  adapted  to  { Ft }.  Let  TZ+  be  the  cr-field  generated  by 
set  [0,  oo).  A  function  /  from  (fi,  Ht)  into  ([0,oo),  TZ+)  is  measurable  if  for  every 

5  €  f-\s)  £  nt. 

Consider  that  neuron  A  is  influenced  by  a  family  of  neurons  B{,  i  =  1, 2,  •  •  • ,  n. 
The  model  we  use  is  depicted  in  Fig.  1;  it  is  similar  in  many  respects  to  that 
studied  by  Knox  (1974)  and  by  Van  Den  Boogaared  et  al.  (1986).  A  sequence  of 
impulses  from  neuron  B{  is  transformed  into  a  membrane  potential  in  neuron  A.  If 
the  integrated  membrane  potential  exceeds  a  threshold  value  6{t),  an  impulse  (spike) 
is  generated  while  the  membrane  potential  discharges  to  a  resting  level  va.  hi{t,s) 
is  the  impulse  response  (not  necessarily  time-invariant)  which  describes  the  total 
temporal  influence  of  neuron  Bi  on  neuron  A  from  past  up  to  present,  including  the 
conduction  and  transmission  delay.  A  synaptic  connection  is  said  to  be  excitatory  if 
h(t,  s)  >  0  for  all  t,  all  s  in  R;  it  is  said  to  be  inhibitory  if  h(t,  s)  <  0  for  all  t,  all  s 
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in  R. 


The  somatic  (membrane)  potential  (WA)  of  neuron  A  is  represented  by  a  linear 
spatial-temporal  superposition  of  all  input  action  potentials  of  neurons  B\ ,  i?2>  •  •  • ,  Bn 
(including  self-inhibition  and/or  self-excitation),  and  an  unknown  random  potential 
Ut  which  represents  the  influence  of  all  other  unobservable  neurons  and  biophysical 
factors.  A  sigmoid  function  g  is  used  to  map  the  somatic  potential  as  follows: 

t 

WA  =  g(Ut  +  £  /  h(t,r)dNBi(r))  =  g(Ut  +  £  £  k(t,T^))  (1) 

i  J0  i  k= 1 

where  {T^'}k>  1  are  the  epoch  times  of  spike  train  from  neuron  B{,  and  {NBi(t)}t> o 
is  the  associated  counting  process,  i.e.,  the  number  of  spikes  arriving  from  neuron  /?, 
in  the  interval  (0,  <]- 

For  mathematical  simplicity,  let  us  assume  that  the  nonlinearity  g  has  the  form  of 
(7(2:)  =  aex,  a  >  0,  i.e.,  that  neuron  A  is  operating  around  threshold  and  is  thus  not 
strongly  driven.  Suppose  further,  without  loss  of  generality,  that  we  are  interested 
in  finding  the  connectivity  between  two  neurons  ( A  and  Bi).  Then  we  write 

wtA  =  -  g(ut  +  £  £  h(t,T*))  vtA  (2) 

»/ 1  fc=i 

where  VA  is  called  here  the  semi-membrane  potential  and  is  defined  as: 

NBl  (■ t ) 

VtA  =  g(  £  Kt,T*')).  (3) 

fc= 1 

In  order  to  account  for  the  firings  of  neuron  A  that  are  due  to  VtA ,  we  can  think  of 
the  factor  %g(Ut  +  h(t,T^1))  as  a  continuous  random  variable  Z,  such 
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that  a  random  threshold  9(t)  is  formed,  which  is  defined  as: 


o(t)  =  ze0{t)  (4) 

where  90(t)  is  a  time  function.  Due  to  the  refractory  period  r  during  which  a  neuron 
is  unable  to  produce  a  successive  spike,  the  time  function  can  be  taken  as  simple  as 

f  oo,  Tk  <  t  <  Tk  +  r 

«.(<)  =  {  (5) 

(  Oo,  Tk  +  r  <  t  <  Tk+i 

where  90  >  0  is  a  constant,  Tk  and  Tk+ 1  are  the  times  at  which  the  fc-th  and  (k  +  l)-st 
spike  occur,  respectively. 

A  spike  occurs  whenever  the  threshold  is  exceeded  by  the  accumulated  semi¬ 
membrane  potential,  i.e., 

fvTAdr>e{t)  (6) 

Jto 

where  to  is  the  instant  of  the  preceding  spike.  Denoting  by  iVyi(t)  the  number  of 
spikes  in  train  A  during  time  interval  (0,  /],  a  stochastic  counting  process  {iV^(t)}f>0 
is  associated  with  spike  train  A  with  Na( 0)  =  0.  Let  AiV^t)  =  NA(t  -f  At)  —  N^^t) 
be  the  number  of  spikes  in  an  infinitesimal  duration  At.  We  say  that  a  process  is 
orderly  if  Pt(ANa(T )  >  1)  =  o(At). 

Armed  with  this  general  neuron  model,  we  are  ready  for  the  analysis  of  the 
interneuronal  connectivities  deduced  from  the  stochastic  firing  of  several  neurons. 


3  Analytical  Results 

In  this  section,  we  shall  derive  and  elaborate  on  three  basic  results.  We  will  first 
consider  the  simpler  case  of  two  observable  neurons,  and  show  how  the  connectivity 
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between  them  can  be  expressed  analytically  in  terms  of  the  neuron  model  outlined 
above.  We  then  consider  the  sources  of  uncertainty  in  this  estimate  and  how  they  can 
be  reduced  through  added  information  from  neighboring  neurons.  Finally,  we  will 
comment  on  the  critical  normalization  procedures  used  to  remove  the  confounding 
effects  of  stimulus  artifacts.  In  the  following  discussion,  assume  that  At  =  {  neuron 
A  firing  at  instant  i},  and  Bs  =  {  neuron  B  firing  at  instant  5}.  The  three  basic 
results  derived  are: 


Result  1.  The  joint  firing  probability  of  a  pair  of  presynaptic  and  post-synaptic 
neurons  can  be  expressed  as  the  product  of  individual  firing  probabilities  and  the 
pairwise  connectivity,  and  a  corrupting  (uncertainty)  factor  due  to  other  unobserv¬ 
able  influences  on  the  firing  of  A: 


Pr{AuBa)  =  Pr(At)Pr(Bs)1(t,s)eh 

where  7 (f,  s)  is  the  corrupting  factor  (7  >  0)  given  by: 

.  _  E[fA(tA)PB(s)\ 
E[fA(t,9t))E[PB(VB)Y 


(7) 


(8) 


where  ./U(t,#t)  is  measurable  with  respect  to  , Na(t);t  <  t )  and  Pb(s )  is 

measurable  with  respect  to  a(V^,  Nb(t);  t  <  s )  . 


Result  2.  The  uncertainty  can  be  reduced  (i.e.,  the  corrupting  factor  made  closer  to 
1)  if  more  interacting  neurons  are  observed  simultaneously.  The  pairwise  connectivity 
becomes: 


h(t,  s)  =  log 


Pr(At,Bs,Ct)Pr(Ct) 

Pr(At,Ct)Pr(Bs,Ct ) 


-log  7* 


(9) 
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with  Ct  =  {neurons  C\,  C2,  ■  •  -,Cm  firing  at  instant  f},  where 


~(i  ,  =  E[fA(t,Ot)PB(s)\Ct} 

7  U  J  E[fA(t,dt)\Ct]E[PB(s)\Ct] 


(10) 


is  a  quantity  satisfying  |q*  —  1|  <  |q  -  1|.  If  7*  is  very  close  to  1,  then  log  7*  may  be 
negligible. 


Result  3.  In  order  to  minimize  the  effects  of  the  stimulus  on  the  estimators  of  the 
connectivity,  the  normalized  joint  firing  probability  given  by  : 

N„(t,  s)  =  Pr(At,Bs)/Pr(At)Pr(Bs)  (11) 

leads  to  estimators  superior  to  those  produced  by  the  often  employed  shuffle  method: 

N3(t,s)  =  Pr(At,Bs)  -  Pr(At)Pr(Bs),  (12) 

which  is  the  quantity  that  the  cross-covariance  histogram  estimates. 

3.1  Definitions  and  Further  Relationships 


In  order  to  discuss  the  derivation  of  the  above  stated  results,  we  will  need  to 
utilize  a  few  more  definitions  and  relationships.  In  particular,  we  will  make  use  of 
the  PST  histogram  of  a  single  cell  spike  train  which  measures  the  firing  rate  of  a 
neuron  with  respect  to  the  stimulus  onset.  Each  bin  of  the  PST  histogram  is  an 
unbiased  estimator  for  the  probability  density  of  the  average  neuron  firing  over  a 
short  period  At  at  instant  t  corresponding  to  that  bin.  The  firing  probability  density 
of  a  neuron  (A)  is  defined  as 

v.m  ■ 

where  Pa( t)  is  the  conditional  firing  probability  density  of  the  post-synaptic  neuron 
given  the  history  of  the  intensity  process  of  the  presynaptic  neuron  and  the  history 
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of  spike  train  A  denoted  by  MA  —  a(NA(s)  :  s  <  t),  that  is, 


=  lim  ■  m s-m. 

V  '  At— *0  At 

Therefore,  the  PST  histogram  of  a  neuron  A  estimates  E[PA(t)]At. 
Likewise,  we  have 


(14) 


-  j&J,r(AJA-,)  =  1)  <») 

where  Pb(&)  is  the  conditional  firing  probability  density  of  the  presynaptic  neuron 
given  the  history  of  the  intensity  process  of  the  presynaptic  neuron  and  the  history 
of  spike  train  B,  that  is, 


Pb(s) 


lim  Pr(ANB(s)  =  l\nf,AfsB) 
As— >0  As 


(16) 


Therefore,  the  individual  PST  histograms  of  neurons  A  and  B  estimate  E[PA(t)]At 
and  J5[Pb(s)]As,  respectively. 


Moreover,  the  joint  PST  histogram  of  the  two  neurons  estimates  E[Pab{^i  s)]AtAs 
where 


E[PAB(t,s)\ 


lim  Pr(^NA(t)  =  l,ANB(s)=:l) 
A<,A«-K)  At  As 


(17) 


and 


PAB{t,s ) 


lim 

At, As— f0 


Pr(ANA(t)  =  1,  AJVS(«) 


At  As 


(18) 


Given  a  pair  of  interacting  neurons  ( A  and  B),  the  following  lemmas  will  play  an 
important  role  in  the  analysis  below. 
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Lemma  1.  PB (f)  can  be  expressed  as  a  map  from  the  semi-membrane  potential  space 
of  neuron  B  onto  [0,  oo), 

/bOQ  vb 

1  -  FB(xt)  1 

where  xt  =  j£  Vfdr,  and  fB(-),  FB{-)  are  the  density  and  the  distribution  functions 
of  the  threshold  of  neuron  B,  respectively.  The  function  PB(-)  may  have  a  very  simple 
form.  For  example,  if  the  threshold  is  an  exponentially  distributed  independent 
random  variable  with  mean  A,  then  PB(t)  =  AFtB .  And  in  this  case,  {NB(t)}t> o  is  a 
doubly  stochastic  Poisson  process. 

Lemma  2.  The  conditional  expectation  of  the  product  of  the  semi-membrane  poten¬ 
tial  of  neuron  A  and  the  firing  rate  of  neuron  B  can  be  expressed  as 

E[vAdN^)}  =  eh(t,s)E[yApB{s)]  (20) 

The  proofs  of  the  lemmas  are  given  in  the  Appendix.  Lemma  1  gives  the  ex¬ 
pression  of  the  conditional  firing  probability  density  of  the  presynaptic  neuron  given 
the  history  of  the  intensity  process  of  that  neuron.  Lemma  2  relates  the  connectiv¬ 
ity,  the  membrane  potential  of  the  post-synaptic  neuron,  and  the  firing  rate  of  the 
presynaptic  neuron. 

3.2  Discussion  of  Result  1 

We  will  first  need  to  derive  an  expression  for  the  firing  rate  of  the  post-synaptic 
neuron  (A).  In  general,  the  threshold  9t  of  this  neuron  is  not  an  independent  variable, 
since  it  depends  on  all  other  unobservable  inputs  to  the  neuron.  Given  an  arbitrary 
value  for  6t,  we  can  write 

Pr(ANA(t)  =  l|0t)  ~ 
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(21) 


rt+At  a 

Pr( l  VTAdr  >  0t\0t)  ~  Pr(VtA  >  jL\9t). 

Since  6t  is  assumed  to  be  a  positive  threshold,  by  the  Markov  inequality  we  have 

Pr(ANA(t)  =  l|0t)  <  ^E[VA\Ot\.  (22) 

Vt 

Averaging  for  all  possible  6t  and  taking  limit  as  At  goes  to  zero  result 

Pr{At)  =  Pr{dNA(t )  =  1)  <  E[~VA]dt.  (23) 

In  fact,  it  can  be  shown  (see  the  Appendix)  that  there  exists  a  measurable  function 
fA(t, 6t)  with  respect  to  a (Va,Na{t)\t  <  t)  such  that 

Pr(At)  =  E[fA(t,6t)]dt.  (24) 

Similarly,  the  conditional  joint  firing  probability  can  be  expressed  as 
Pr(ANA(t)  =  1,ANb(s)  =  l\9t)  a 

rt+At  At  , 

PA  jf  VA dr  >  9U  ANb(s)  =  l|0t)  <  ~E[VAANB(s)\0t).  (25) 

By  Lemmas  1  and  2,  we  therefore  have 

Pr(Au  Bs )  =  Pr{dNA{t)  =  l,dNe(s)  =  1)  <  eh^ E[^VAPB(s)]dtds.  (26) 

tq 

As  in  Eq.(24)  above,  a  more  precise  expression  can  be  derived  as 

Pr(At,  B.)  =  eh^E[fA(t,  et)PB{s))dtds.  (27) 

Since  the  firing  probability  of  the  presynaptic  neuron  is 

Pr(B. )  =  E[PB(s))ds,  (28) 
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then  combining  equations  (27),  (24)  and  (28)  gives  Result  1  with 


_  E[fA(t,Ot)PB(s)} 
E[fA{tA )]  E[Pb(s)Y 


(29) 


3.3  Discussion  of  Result  2 

The  factor  7 (t,s)  reflects  our  ignorance  of  the  input  to  neuron  B,  or  that  of 
the  knowledge  of  the  threshold  0t.  For  a  completely  known  input  {V^8}  (hence 
Pg(s)  is  determined),  7 (t,s)  —  1;  for  a  completely  known  threshold,  7 (t,s)  is  a 
constant.  When  the  activity  of  more  neurons  are  known,  the  uncertainty  in  the  input 
and/or  threshold  decreases,  and  j(t,s)  approaches  1.  For  instance,  if  the  activities 
of  more  interacting  neurons  (Ci,  C2,  •  •  • ,  Cm)  are  available,  we  can  use  a  multiunit 
PST  histogram  in  addition  to  the  conventional  joint  and  individual  histograms  to 
estimate 

Pr(At,  Bs,  Ct)F‘r(Ct)  _  PrjAt,  Bs\Ct)  _  \  h(t,s )  /  oa\ 

Pr(At,Ct)Pr(Ba,Ct)  Pr(At\Ct)Pr(Bs\Ct)  7  K  '  ^ 

where  7 is  defined  in  Eq.(10). 

Note  that  observing  more  neurons  makes  7*  closer  to  one  than  7  is  in  Eq.(29), 
and  hence  the  estimator  for  h(t,  s )  is  more  reliable.  This  is  because  neurons  Ct  may 
contain  information  about  Pb(s)  and/or  6t  (for  instance,  if  neurons  Ct  influence  the 
activity  of  either  or  both  neurons  A  and  B),  hence  making  /^(t,^)  less  correlated 
with  Pe(s),  and  7*  closer  to  one. 

3.4  Discussion  of  Result  3 

An  important  factor  in  correctly  interpreting  the  correlations  among  the  activities 
of  different  cells  concerns  the  effects  of  the  stimulus.  Specifically,  this  refers  to  the 
fact  that  unconnected  cells  may  exhibit  strong  correlations  in  their  firings  purely 
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due  to  the  fact  that  they  are  driven  by  the  same  stimulus.  In  order  to  eliminate 
these  effects,  some  form  of  normalization  is  necessary.  In  Result  3  we  show  how 
the  stimulus  shuffle  alone  fails  to  accomplish  this  task.  In  order  to  illustrate  this 
with  explicit  analytic  expressions,  three  simplifying  assumptions  will  be  adopted 
concerning  the  properties  of  the  post-synaptic  neuron  threshold  6t  (used  in  Theorem 
1  below)  and  the  distribution  of  the  presynaptic  potential  (used  in  Theorem  3).  We 
start  by  stating  two  of  these  assumptions  and  the  theorems  associated  with  them, 
and  then  proceed  to  relate  the  correlation  functions  explicitly  to  the  inter-neuronal 
connectivity  (h(t,  s))  in  a  pair  of  neurons  ( A  and  B ). 


Assumption  1.  The  random  variable  Z  of  the  threshold  in  Eq.(4)  is  independent  of 
VA,  and  has  an  exponential  pdf: 


Pz(*) 


0oe~(9oZ~Vo\  vo/0o  <  z  <  oo 

< 

k  0,  z  <  v0/90  . 


(31) 


This  assumption  is  typically  valid  in  cases  where  Neuron  B  is  only  related  to 
neuron  A,  i.e.,  it  is  weakly  related  to  any  other  neuron.  It  can  be  verified  that  under 
Assumption  1,  the  output  spike  train  of  the  neuron  A  is  a  doubly  stochastic  Poisson 
process  {Ar/t(t)}j>o  with  the  intensity  process  (Aj*}j>o,  where 


A' 


0,  Tk  <  t  <  Tk  +  r 
.  VtA,  Tk  +  r  <t  <Tk+ 1 


where  r  is  the  refractory  period. 


(32) 


For  a  doubly  stochastic  Poisson  process  {^^(t)},  we  have 

Pr(ANA(t)  =  0\H?,AftA)  =  l-AAAt  +  o(At), 
Pr(ANA(t)  =  l\n^KA)  =  A?At  +  o(At), 
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(33) 

(34) 


Pr(ANA(t)>l\K?,AftA)  =  o(At)  (35) 

Note  that  depends  on  A ftA,  and  hence  {NA(t)}t>o  is  a  self-exciting  process  with 
the  intensity  function  E[Af\AfA]. 

Assumption  2.  The  refractory  period  is  much  smaller  than  any  interspike  interval 
and  hence  is  negligible. 

Under  this  assumption  Af  does  not  depend  on  AftA,  and  hence  the  intensity 
process  becomes  the  membrane  potential  process. 

Theorem  1.  Under  Assumptions  1  and  2,  {NA(t)}t>o  is  a  doubly  stochastic  Poisson 
process  with  the  intensity  process  {VA}t>o-  Furthermore,  the  conditional  joint  firing 
probability  density  of  neurons  A  and  B  can  be  expressed  as 

PAB(t,s)  =  PA(t)PB(s)eh ■<*••)  (36) 

for  all  t  and  all  s,  where  h(t,s )  is  the  inter-neuronal  connectivity  with  non-zero 
transmission  delay. 

Theorem  1  states  that  the  joint  firing  probability  density  can  be  expressed  as 
the  product  of  the  individual  firing  densities  and  the  connectivity.  Thus  the  inter¬ 
neuronal  connectivity  h(t,  s)  can  be  directly  identified  by 

h(t,s)  =  log PAB(t,s)  -log PA(t)  -log PB(s).  (37) 

This  is  an  ideal  case. 

Experimently,  if  only  neurons  A  and  B  are  recorded,  the  semi-membrane  po¬ 
tential  of  the  presynaptic  neuron  B  is  generally  unknown  (hence  PAB ,  PA,  and  PB 
are  unknown),  because  the  membrane  potentials  are  unobservable  in  extracellular 
recordings.  Therefore,  one  must  instead  evaluate  the  normalized  unconditional  joint 
probability  density  Np(t,s). 
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Let  us  recall  that  the  PST  histogram  of  neuron  A  estimates  E[PA(t)\At  and 
that  of  neuron  B  estimates  E[Ps(t)]At,  and  the  joint  PST  histogram  estimates 
E[PAB(t,  s)\At As.  Therefore,  Np(t,s )  can  be  formed  by  these  three  histograms: 


Np(t,s ) 


E[PAB(t,s)] 

E[PA(t)]E[PB(s)Y 


(38) 


Lemma  3.  Under  Assumptions  1  and  2,  the  conditional  firing  probability  of  the 
post-synaptic  neuron  can  be  expressed  as 


PA(t)  =  E{VtA\Hf]  =  a  exp 


h(t,r) 


1  )V*dT}. 


(39) 


Theorem  2.  If  Assumptions  1  and  2  hold,  the  uncertainty  7  in  Result  1  can  be 
expressed  as 


E[VB]  fJ[exp{J2(cM*.r)  -  l)VBdr}]' 


(40) 


Theorem  3.  If  the  semi-membrane  potential  process  of  the  presynaptic  neuron  can 
be  decomposed  in  the  form  of  VtB  =  Xf(t )  where  X  is  a  random  variable,  and  /(f) 
is  a  deterministic  time  function,  then 


7(M)  = 


M'(yt) 
E[X]  M{j]t) 


(41) 


where  M(-)  is  denoted  as  the  moment  generating  function  of  X,  and  M'(rjt )  is  the 
first  derivative  of  M(-)  with  respect  to  rjt  which  is  expressed  by 

m  =  / V(t,T)  -  1  )f(r)dr.  (42) 

Jo 

Furthermore,  i(t,5)  — >  1  when  V ar(X)  — >  0. 


16 


Assumption  3.  The  presynaptic  membrane  potential  X  is  Gamma  distributed  with 
parameters  (A,i/). 

One  consequence  of  Theorem  3  is  that  if  Assumption  3  holds  —  a  relatively 
common  occurrence  (Bishop  et  al.  1964;  Nakahama  et  al.  1968;  Correia  and  Lan- 
dolt  1977)  —  the  normalized  unconditional  joint  probability  density  Np(t,s )  can  be 
explicitly  evaluated  in  terms  of  these  parameters  as 

Np(t,s)  —  — eh(t,a\  m  <  A.  (43) 

Comparing  this  expression  with  that  of  Theorem  2  suggests  that  7 (t,  s )  =  A/(A  —  rjf). 
Therefore,  for  a  given  Gamma  distribution  (of  degree  u),  as  the  variance  of  X  (— 
v/X2)  becomes  smaller,  A  increases,  and  7 (t,s)  —>  1.  In  other  words,  the  more  is 
known  about  VSB  (e.g.,  from  recordings  of  additional  neurons),  the  more  accurate  is 
the  estimate  of  the  connectivity  between  neurons  A  and  B.  We  will  illustrate  these 
result  through  simulations  later  in  section  5  (see  Fig.  3). 

If  the  membrane  potential  does  not  vary  much  for  different  stimulus  presentation 
(small  variance  of  X),  then  A  >  fft.  Consequently,  we  have 

Np(t,s)  ~  eh(t,a\  (44) 

This  confirms  the  conclusions  established  in  Result  2  earlier. 

In  contrast  to  the  normalization  used  in  Eq.(ll),  the  conventional  cross-covariance 
histogram  (which  is  the  modified  joint  PST  diagram  using  the  shuffling  method)  uses 
a  difference  normalization  which  estimates  (Gerstein  1970;  Habib  and  Sen  1985) 

Nd(t,  s )  =  E[PAB(t,  ^)]  -  E[PA(t)]  E{Pb(s)}.  (45) 

In  general,  this  expression  is  very  complicated.  However,  if  we  make  use  of  the 
assumptions  in  Theorem  3  for  the  intensity  process  of  the  presynaptic  neuron  (i.e., 
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a  Gamma  distribution),  it  reduces  to 


Nd(t,s )  =  av  f  (s)M  (r)t)(- - -).  (46) 

A  —  T]t  A 

If  the  membrane  potential  is  not  varying  too  much  for  different  stimulus  presentation 
(A  rjt),  then  N<t(t,s)  can  be  approximately  written  as 

Nd(t,  s )  ~  (47) 

This  expression  suggests  that  identifying  the  connectivity  here  is  considerably  more 
difficult  than  that  of  the  normalization  Np(t,s )  used  earlier,  since  quantities  a,  n,  A 
and  function  f(s )  are  generally  unknown.  Nevertheless,  Eq.(47)  suggests  that  the 
shuffling  method  remains  effective  in  indicating  the  absence  of  a  direct  connection 
(i.e.,  when  h(t,s)  is  very  small),  since  in  that  case  Nd(t,s )  is  approximately  zero 
regardless  of  the  confounding  terms  (a,  v,  A  and  function  /(s)).  We  will  illustrate 
this  in  simulations  in  section  5  (see  Figs.  4). 


4  Experimental  Considerations 

In  the  analysis  of  multi-neuronal  connectivities,  spike  trains  from  several  neurons 
are  recorded  in  response  to  the  repeated  presentation  (e.g.,  R  times)  of  a  stimulus. 
Spikes  are  usually  sampled  and  parsed  into  (i.e.,  labeled  by)  small  time  bins,  using 
the  onset  of  the  stimulus  as  the  initial  bin.  The  bin  width  At  is  always  chosen 
to  be  so  small  that  at  most  one  spike  may  occur  in  each  bin  (which  corresponds 
to  the  orderliness  of  the  point  process).  Thus  each  spike  train  is  converted  into  a 
discrete  0-1  process,  and  is  further  segmented  into  R  segments,  each  for  one  stimulus 
presentation. 
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Let  Arn  be  the  time  bin  corresponding  to  the  n-th  bin  associated  with  the  r-th 
stimulus  presentation.  A  spike  train  can  then  be  represented  by  a  R  X  N  random 
matrix  A  with  elements  (Arn,  r  —  1,2,***,  jR;  n  =  1, 2,  •  •  ■ ,  N)  —  called  here  a  spike 
matrix.  Let  us  assume  that  the  firing  activity  during  each  stimulus  presentation 
is  statistically  independent.  Therefore,  each  element  is  a  random  variable  taking 
values  {0,1},  and  the  elements  in  the  same  column  are  independent  and  identically 
distributed. 

The  PST  histogram  ( Hn ,  n  =  1, 2,  •  •  • ,  N)  reflects  the  stimulus-locked  firing  rate 
of  each  single  neuron,  and  it  is  formed  by  taking  average  over  every  column  of  the 
spike  matrix, 

1  R 

Hn  =  n  =  1,2,  ••  -,1V.  (48) 

r= 1 

The  value  of  H£  counts  the  average  spikes  over  R  stimulus  presentations  in  the  n-th 
bin  in  a  spike  train  A. 

The  joint  PST  scatter  diagram  of  two  neurons  A  and  B  m  =  1,2  n  — 

1, 2,  •  •  • ,  N)  measures  the  coincidence  spikes  in  train  A  and  in  train  B  relative  to  stim¬ 
ulus  onset.  It  is  a  two-dimensional  histogram  with  one  axis  (m)  for  train  A  and  the 
other  axis  (n)  for  train  B,  and  hence  it  is  an  N  square  matrix  H.  Element  H 
represents  the  average  count  for  coincidence  of  a  spike  in  the  m-th  bin  of  train  A  and 
a  spike  in  the  n-th  bin  of  train  B  over  R  stimulus  presentations,  that  is, 

1  R 

Hmn  =  ji'EArmBrn,  m  =  1.,  2  *  *  * ,  JV ;  n  =  1.,  2,  *  *  - ,  iV  (49) 

r=l 

where  Arm  and  Brn  are  the  elements  of  spike  matrices  for  trains  A  and  B,  respectively. 
Therefore,  the  matrix  presentation  of  the  joint  PST  scatter  diagram  is 

H  =  4aTb  (50) 

K 
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where  T  denotes  transposition.  The  expanded  joint  PST  histogram  for  multiunit 
recordings  (of  M  neurons)  is  then 

1  R 

lining  =  ^C1rniC^2---C^M,  nt  =  l,2-..,iV;  t=l,2,--.,M  (51) 

L  r= 1 

where  Clrn.  is  the  element  of  the  spike  matrix  for  the  i-th  neuron. 


4.1  Using  the  Scatter  Plot  to  determine  neuronal  connectivities 


The  correlations  between  a  pair  of  recorded  neurons  (A  and  B )  can  be  computed 
from  the  experimental  estimate  of  the  expression  of  Result  1,  i.e., 


/>((,*)  =  log(^ltil)  =  log( 


E[PAB(.t,s)\ 
E[PA(t)}  E[PB{s )] 


)  ~  l°g(7(i,  ■«)) 


where  and  E[Pb(s)]  represent  the  PST  histograms  of  firings  of  the  neuron 

pair,  E[PAB(t,s)]  is  their  scatter  plot,  and  7  (>  0)  is  the  corrupting  factor  represent¬ 
ing  the  uncertainty  in  the  estimate  due  to  the  influences  of  other  unobserved  neurons 
and  biophysical  factors.  Thus  in  terms  of  bin  numbers  m  and  n,  the  above  equation 
can  be  written  as 


Hab 

l°g =  h(m&t,  nAt)  +  \og(-y(mAt,  nAt)).  (52) 


In  the  case  of  time  invariant  connectivities,  h(t,s )  becomes  h{t  —  s),  and  the  corre¬ 
lation  peak  becomes  a  band  that  runs  parallel  to  the  principal  diagonal  (t  —  s  =  0).1 


In  the  practical  application  of  Eq.(52),  the  confounding  7 (t,s)  contributions  are 
not  known.  However,  the  analysis  shows  that  additional  simultaneous  recordings  can 


1Note  that  one  can  detect  further  correlations  in  the  unnormalized  scatter  plot,  such  as  the 
more  diffuse  bands  of  time-invariant  common  inputs  (Yang  1989).  Of  course,  these  features  are 
intentionally  removed  by  the  normalization  since  they  do  not  reflect  direct  connectivities  between 
the  neuron  pair. 
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be  used  to  reduce  these  uncertainties.  Therefore,  by  using  the  additional  data,  the 
improved  estimator  for  h(t,  5)  becomes 

jjABC3-Cm  jjC3-Cm 

log(^^TF BflggJ  =  M™Af,nAf)  +  log(7*(mAf,nAt))  (53) 

where  are  s^mPly  the  joint  multi-dimensional  scatter  plots  defined  in 

Eq.(51),  and  the  uncertainty  factor  7*  (<  7)  is  defined  in  Eq.  (10).  The  estimates 
of  Eqs.(52)  and  (53)  are  illustrated  in  network  simulations  in  section  5. 

4.2  Establishing  confidence  measures  on  the  estimates 

The  histograms  are  random  variables  subject  to  fluctuations.  Hence,  it  is  impor¬ 
tant  to  determine  upper  and  lower  bounds  such  that  we  assume  a  connection  between 
neurons  A  and  B  whenever  these  bounds  are  surpassed.  By  the  law  of  large  numbers, 
converges  to  E[PAB(t,s)\,  so  does  to  E[PA(t )]  and  11%  to  E[Pb(s)]  almost 
surely  as  R  — >  00.  Therefore,  if  neurons  A  and  B  are  independent,  by  theorems  on 
limiting  distributions, 

jf'FfjB  1  as  R  00  (54) 

11  m  n 

almost  surely. 

The  hypothesis  Bo  is  that  the  two  neurons  are  statistically  independent,  which 
is  supported  by 

E[PAB(t,  5)]  =  E[PA(t))  E[Pb(s)}.  (55) 

And  the  alternative  hypothesis  B\  is  that  the  two  neurons  depend,  which  is  described 

by 

E[PAB(t,  5)]  ^  E[PA(t)}  E[Pb(s)}.  (56) 

One  expects  /  H%II%  to  be  close  to  1  if  hypothesis  Bo  is  true.  Conversely,  if 
the  amount  it  deviates  from  1  exceeds  a  bound  b,  one  accepts  hypothesis  B\. 
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Now  for  a  given  significance  level  a ,  we  need  to  find  the  bound  b  satisfying 


Pr(  I 


hab 

HAHB 


-l\>b\n0)  =  a. 


The  hypothesis  testing  is  stated  as  the  following  theorem. 


(57) 


Theorem  4.  Let  b  be  a  bound  which  divides  a  critical  region  for  the  hypothesis 
testing.  One  announces  that  there  is  a  dependence  between  the  two  observed  neurons 
if 


H 


AB 


HAHB 


1|  >  b. 


(58) 


For  the  given  significance  level  a  of  false  announcement  of  dependence,  the  bound 
can  be  calculated  by 


b  =  £b 


i  -b&h* 

RHAHB 


(59) 


where  the  value  of  £(,  is  determined  from 


*(66)  =  1  ~  | 


(60) 


and  $(®)  =  ^S-ooe  X?,2dx. 


The  function  $(a;)  is  usually  available  as  the  standard  normal  distribution  table. 
For  example,  a  =  0.05  gives  £b  =  1.96. 


The  above  theorem  implies  that  element  HAB/(HAHB)  of  the  normalized  joint 
PST  diagram  has  a  conditional  expectation  value  1  and  an  approximate  conditional 
variance 


a 


2 

mn 


rv»/ 


1  -  Ha Hb 
RHA II B 


(61) 
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given  the  values  of  II £  and  under  hypothesis  Hq.  Since  H*  and  H, f  are  usually 
very  small  and  R  is  fairly  large,  this  approximation  is  close  to  a  recent  result  by  Palm 
et  al.  (1988)  where  their  conditional  variance  is 


.2  (1  ~  Hj){  1  -  Hi) 

mn  (R-1)H*H2 


(62) 


under  hypothesis  Ho- 

The  bound  dividing  the  hypothesis  regions  can  be  made  more  useful  in  neural 
networks  with  time  invariant  connectivities.  Let  wm  reflect  the  fluctuation  in  the 
normalized  joint  PST  diagrams  such  that 

JfAB 

jTjfs  =  7(™Af ,  nAt)eh«m~' ^  +  wm  ,  (63) 

and  the  mean  of  wn  is  zero.  Let  k  =  m  —  n.  A  collapsed  version  can  be  generated 
by  averaging  over  diagonals  of  the  normalized  joint  PST  diagram.  This  collapsed 
version  is  a  1-dimensional  histogram  Gk  expressed  by 

1  min  (N,N-k)  fjAB 

s-i  _  v-'  Jin+fc,n 

‘  ”  N  -  W  n=Jw-*>  H^S 

k  =  -N  +  1,  0,1,  ■■■,!!-  1  (64) 

where  k  =  0  is  the  collapsed  point  of  the  principal  diagonal. 


Since  averaging  reduces  the  fluctuations  (the  average  of  wm  has  a  smaller  vari¬ 
ance),  Gk  is  a  better  estimator  for  the  time  invariant  connectivity  h(t,s)  =  h(t  —  s ). 
This  enables  us  to  establish  a  bound  such  that 

PT(\Gk  -  1|  >  bk  |  Ho)  =  a.  (65) 
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Theorem  5.  Given  a  significance  level  a,  let  b be  a  bound  of  critical  region  satis¬ 
fying  the  above  equation,  then  b k  may  be  approximately  written  as 


h  - 


min(N  ,N —k) 
n=max(l,l-fc) 


a 


2 

n+k,n 


N  -  \k\ 


£6 


(66) 


where  £{,  is  the  the  same  as  in  Theorem  4,  and  o^n  is  given  by  (61).  Furthermore, 
bk  will  reduce  to 

.  b 
k  JN^W\ 

when  o'*'  s  are  taken  as  constants. 


This  theorem  indicates  that  the  critical  region  is  enlarged  (the  bound  value  de¬ 
creases)  when  the  collapsed  version  of  the  normalized  joint  PST  histogram  is  used. 


5  Simulations  and  Discussion 

In  order  to  illustrate  the  nature  of  the  estimates,  uncertainties,  and  bounds  de¬ 
rived  earlier,  we  show  the  results  from  simulations  of  networks  of  excitatory  and 
inhibitory  neurons.  The  neuron  model  used  for  the  simulations  is  depicted  in  Fig. 
1(c)  where  the  nonlinearity  g(x )  =  aex  and  the  random  threshold  has  an  exponential 
distribution  with  mean  1. 

In  the  first  case  (Fig.  2),  pair-wise  excitatory  and  inhibitory,  time-invariant 
connections  are  estimated  using  the  normalized  scatter  plots;  the  uncertainty  factor 
(7)  is  equal  to  1.  The  upper  plots  show  the  two-dimensional  normalized  scatter 
plots.  The  correlations  appear  as  bands  along  the  principal  diagonal  because  h(t,  s) 
is  time-invariant.  Hence,  the  scatter  plot  can  be  collapsed  along  this  axis  to  produce 
the  lower  histograms.  Note  that  time- variations  in  h(t,s)  (e.g.,  due  to  post-stimulus 
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adaptation)  do  not  allow  this  reduction.  Consequently,  it  should  only  be  performed 
on  the  portions  of  the  neural  record  that  display  obvious  stationary  behavior.  In 
both  simulations  of  Fig.  2,  the  predicted  analytical  estimates  are  also  plotted  for 
comparison,  together  with  the  bound  lines  for  the  confidence  measures  (determined 
by  Theorem  5). 

In  order  to  illustrate  the  effects  of  the  uncertainty  factor  7,  we  examine  in  Fig. 
3  the  interactions  among  three  neurons  with  time-invariant  connectivities.  Here, 
neuron  A  is  inhibited  by  neuron  B  and  excited  by  neuron  C,  and  neuron  C  is  in  turn 
excited  by  neuron  B.  Because  of  the  interactions  between  B  and  C,  the  threshold 
in  neuron  A  is  no  longer  independent  of  the  firings  of  B.  Thus,  if  we  attempt  to 
identify  the  connectivity  between  neurons  A  and  B  from  pair-wise  recordings,  the 
estimates  will  be  contaminated  by  the  7  uncertainty  factor.  The  top  curve  in  Fig. 
3  first  shows  the  “target”  theoretical  connectivity  obtained  from  the  multi-recording 
estimate  given  by  formula  (30)  with  7 *(t,s)  =  1  (i.e.,  ehAB^t~s'>).  If  neuron  C  is 
ignored,  the  pair-wise  estimate  of  etlAB^t~s'>  is  shown  as  the  middle  curve  in  Fig.  3 
(corresponding  to  formula  (52)).  The  correlation  is  so  distorted  that  actual  inhibition 
becomes  false  excitation  because  of  the  strong  excitatory  activity  from  neuron  C.  To 
correct  the  erroneous  correlation,  we  have  to  use  the  information  from  the  third 
neuron.  The  tripartite  correlation  according  to  formula  (53)  is  displayed  at  bottom 
of  Fig.  3,  which  is  much  closer  to  the  analytical  estimate. 

Fig.  4  compares  the  preferred  normalization  with  the  difference  normalization 
(shuffle  method)  under  two  situations.  In  the  absence  of  a  direct  connection,  the 
shuffle  method  provides  accurate  indication  of  the  lack  of  synaptic  inputs  between 
the  two  neurons.  However,  in  the  presence  of  a  direct  connection,  the  shuffle  method 
fails  to  remove  completely  the  stimulus  correlations  as  indicated  by  the  deviation  from 
the  analytical  results.  Instead,  the  normalization  suggested  in  this  paper  performs 
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well  in  both  cases. 


In  conclusion,  the  above  simulations  confirm  the  proposed  theory.  The  neuron 
model  adopted  is  quite  general  because  (1)  the  synaptic  connectivity  h(t,  s )  represents 
a  time- varying  system;  (2)  the  processes  representing  spike  trains  are  not  necessarily 
Poisson  processes,  and  (3)  the  nonlinear  function  g(x)  —  aex  is  an  approximation 
of  aex/l  +  aex  when  aex  <C  1,  meaning  that  the  neuron  is  operating  at  low  firing 
rates.  Moreover,  our  analytical  Results  1  and  2  do  not  dependent  on  any  further 
assumptions.  Although  the  three  simplifying  assumptions  were  made  in  order  to  see 
Result  3  more  clearly,  we  did  not  use  these  assumptions  in  the  simulations  of  Fig.  4. 

The  analysis  presented  in  this  paper  also  points  to  the  following  sobering  conclu¬ 
sion:  For  multiunit  correlation  analysis  to  play  a  useful  role  in  establishing  the  basic 
circuitry  of  the  nervous  system,  new  technologies  have  to  be  developed  for  stable, 
multi-unit  recordings.  These  requirements  stem  from  the  need  for  extended  simulta¬ 
neous  recordings  from  many  cells  in  order  to  construct  adequate  scatter  histograms 
and  to  minimize  inherent  uncertainty  due  to  unobserved  but  related  activities.  Unfor¬ 
tunately,  neither  of  these  requirements  are  easily  met  at  present,  although  extensive 
efforts  towards  this  goal  are  underway  through  the  use  of  silicon-based  microelectrode 
arrays  (Kruger  1983). 
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Appendix 


Proof  of  Lemma  1:  Since  the  threshold  of  neuron  B  is  a  continuous  random  variable, 
let  its  probability  density  function  and  distribution  function  be  fs(x)  and  FB(x), 
respectively.  Let  Xt  =  //  VBdr,  where  to  is  the  occurrence  instant  of  the  previous 
spike.  From  definition  (16)  we  have 

P  (f\  —  Um  Pr(xt  +  At  ^  ®t\xt  <  A'/^) 

lim  Pr(xt<et<xt+&t\nB,MtB)  _  vtBPdt(xt ) 

At-o  At  Pr(xt  <  et\nB ,MB)  1  ~Fet(xty  {W) 

Furthermore,  if  the  threshold  is  exponentially  distributed  with  mean  A,  then 
—  Fs(xt)  =  A,  and  hence  PB(t)  =  XV B .  In  this  case,  PB(t)  does  not 
depend  on  {^/^(t)},  and  {Ns(t)}t> o  evolves  without  aftereffects. 


Proof  of  Lemma  2:  Because  AN  sit)  can  take  values  0  and  1  only,  by  Eq.(3)  we  have 


-  (.u„v„,V'0M,  T»ui AN  t.\- ,.«»  AfB.ntMgWzlMlAg) 
£/[o!6Xp{  /  J  )}\ANb{s)  l)?Tfvs,A's  J 


(69) 


For  t  >  s,  the  conditional  expectation  in  the  above  equation  can  be  written  as 


NB(t) 


E[a exp{  h(t,T?)}\ANB(s)  = 


A— 1 


=  £[aexp{  53  /i(i,Ef)}exp{/i(t,s  + As)}exp{ 

fc=i 


£ 


/c=A^^(s-f-As)+l 


h(t,TB)}\nB,N 


rB i 


(70) 
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which  becomes 


eMM)  E[VtA\nffATaB]  (71) 

as  As  goes  to  0.  Since  -Ps(s)  is  a  measurable  function  with  respect  to  a (7if  X  A/-®), 
hence 

ElVtA^f~]  =  eh^E[E[VtA\n?,AfsB}PB(s)}  =  eh^  E[VAPB(s)}.  (72) 

For  t  <  s,  we  have 

,Vf]  =  £[F, •*!«?, Vf]  £[^^|Wf,Arf] 

=  (T3) 

hence 

=  E[E[vA{nB  ,^B]PB(5)]  =  EI^PbW],  (74) 

Since  h(t,  s )  represents  a  synaptic  connectivity,  which  is  a  causal  system  with  non¬ 
zero  transmission  delay,  h(t,s )  =  0  for  t  <  s.  Thus  Eq.(71)  holds  for  all  t  and 
s. 


Proof  of  Equation  (24):  Replace  neuron  A  for  neuron  B  in  the  proof  of  Lemma 
1.  Here  we  emphasis  the  threshold  of  neuron  A  so  that  fet(-)  and  re¬ 

place  fB( •)  and  Fb(-),  respectively.  Thus  limAt— o T’r(Ai\r/i(t)  =  l\Uf,J^A)/At  = 
VtAPgt(xt)/l—Fet(xt),  where  Xt  =  f)Q  VAdr.  Although  the  joint  distribution  function 
of  9t  and  Xt  is  unknown  in  general,  we  can  define  a  measurable  function 


fA(t,0t ) 


VtAfetM 

1  -  Fgt{ Xt)' 


(75) 


Therefore,  we  have 

Pr(At)  =  E[Pr(dNA(t)  =  1\Ka ,ArA)\  =  E[fA(t,  0t)]dt.  (76) 
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Proof  of  Theorem  1:  Suppose  that  Assumptions  1  and  2  hold,  and  that  the  threshold 
0t  has  an  exponential  distribution  with  mean  1.  By  the  proof  of  equation  (24)  and 
Lemma  1, 

YimQPr(ANA(t)  =  l\HA,AfA)/At  =  fA(t,Ot)  =  V(A,  (77) 

hence  spike  train  NA(t)  is  represented  by  a  doubly  stochastic  Poisson  process.  There¬ 
fore,  by  Eqs.(32)  and  (34)  we  have 

PA{t)  =  E[VtA\H?].  (78) 

Because  Poisson  process  evolves  without  aftereffects,  the  conditional  probability 
given  the  firing  histories  of  neurons  A  and  B  can  be  split  into 

Pr(ANA(t)  =  1,  AJVb(s)  =  mtKs) 

=  Pr(ANA(t)  =  1| Hf)  Pr(ANB(s)  =  (79) 

By  Eqs.(32)  and  (34),  the  first  factor  is 

Pr(ANA(t)  =  1|  H?)  =  VtAA  t  +  o(At).  (80) 

We  write  the  second  factor  as 

Pr(ANB(s)  =  l\Ht,nfys)  =  E[ANB(s)\nt,nfys],  (81) 

and  we  have 

Pr(ANA(t)  =  ltANs(s)  =  mt'HL)  =  E[VtAANB(s)\nt,Hfya]  +  o(At).  (82) 
By  taking  average  over  the  cr-field  TLA,  we  obtain 

PAB(t,s)  =  E[VtA^j^l\Hfyal  (83) 

which  is,  by  the  proof  of  Lemma  2, 

PAB(t,s )  =  eh^E[VA\H?]PB{s)  =  eh^PA{t)PB(s).  (84) 
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Proof  of  Lemma  3:  Under  Assumptions  1  and  2,  PA(t)  becomes  the  conditional 
expectation  of  the  post-synaptic  membrane  potential,  which  can  be  expressed  as 

NB(t) 

E[VtA\ Wf]  =  E[aexp{  J2  KW* )}|Wf]  (85) 

fc=i 

Define  an  event  Dn  —  {NB(t)  =  n},  which  has  a  conditional  Poisson  distribution 


Pr{Dn\Hf) 


Ujyfdr]’  -r;v‘*r 


n: 


(86) 


Because  {JVb(<)}i>o  is  an  inhomogeneous  Poisson  process  with  the  associated  point 
process  {Tk}k>o  for  the  given  realization  of  the  intensity  process  {UtB}<>0,  it  can  be 
shown  that  (Larson  and  Shubert  1979) 


NB(t) 

£[exp{  £  h(t,TB)}\NB(t)  =  n,W?J  =  [ 


A;=l 


ft  eh(t,T)yBdT 

foVfdr 


(87) 


Therefore,  we  have 


NB(t) 

£[exP{  E  K^Tk)}\NB{t)  =  n,7ff] 


fc= l 


oo 

=  E^exP{  E  Kt,Tk))\NB{t)  =  n,Hf}Pr(Dn\Hf) 

n=0 


oo 


=  E 


(88) 


and  the  summation  in  the  above  equation  is  a  series  expression  for  an  exponential 
function.  Thus  rewriting  it  proves  the  lemma. 


Proof  of  Theorem  2:  By  Lemma  1,  Pb{s)  -  VSB;  by  Theorem  1,  fA(t,0t )  =  VtA.  We 
write 


7  (M)  = 


WAv.B] 

E[VtA}E[VPY 


(89) 
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Note  that  E[VtAVsB]  =  E[E[VA\Hf}VB]  =  E[PA(t)VsB]  and  E[VtA]  =  J5[£[V^|Wf]]  = 
E[PA(t)}.  Hence,  applying  Lemma  3  completes  the  proof. 


Proof  of  Theorem  3:  In  the  expression  (40),  the  numerator  of  7 (t,  s )  can  be  written  as 
E[X f(s)ex,nt}\  the  denominator  can  be  written  as  E[X  f(s)]E[eXnt}.  Eq.(41)  is  true 
by  canceling  f(s).  Furthermore,  as  the  variance  of  X  decays  to  zero,  )  — ►  f-iemt 
and  M{rjt)  — >  emt,  consequently  7 (t,s)  -*  1. 


Proof  of  Theorem  f:  For  a  given  significance  level  a,  we  need  to  find  a  bound  b 
satisfying 

Pr(l^^~ll>bl'Ho)  =  a-  (90) 

Let  us  remember  that  RHAB  is  binomially  distributed  with  parameters 
(R,E[PAB{t,s)])  and  that  HAB  —*•  HAHB  almost  surely  under  TLq.  Hence  by  the 
central  limiting  theorem, 


RIIab  -  RHAHB 
y/RIIAHB(  1  -  11*11%) 


iV(0, 1)  as  R  — ►  00 


(91) 


where  IV(0, 1)  is  denoted  as  a  standard  Gaussian  random  variable.  This  means  that 
Eq.(56)  can  be  approximately  written  as 

Pr(|tf(0, 1)|  >  £j  |  Wo)  =  a  (92) 


where 

bRHjll I 

yjRIIAHB(  1  -  II&HB) 
which  results  in  an  expression  of  the  bound  as 


(93) 


b  -  £b 


1  ~  R-m  II  n 

RB&hb  • 


(94) 
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The  value  of  £f,  is  determined  by 


*(*)  =  1  -  f 


(95) 


where  $(x)  =  ^  e~x^2dx. 


The  above  arguments  imply  that  element  HAB /(II A II B)  of  the  normalized  joint 
PST  diagram  has  a  conditional  expectation  value  1  and  an  approximate  conditional 
variance 


a 


2 

mn 


1  -  II A II B 
RHAHB 


(96) 


under  hypothesis  Ho- 


Proof  of  Theorem  5:  Let  us  note  that  under  hypothesis  Ho,  HAB / IIAHB  is  approx¬ 
imately  Gaussian  distributed  with  mean  1  and  variance  a^n.  Hence 

^  min  (N,N—k) 

I Gk  ~  1|  —  |  ,,  _  j  X^  &n+k,n  Nn(0,  1)|  (97) 

'  n=max(l,l  —  k) 


where  each  1V„(0, 1)  approximately  has  a  standard  Gaussian  distribution  expressed 
by 


Nn(0, 1)  = 


-  unions 


(98) 


and  is  given  in  Eq.(61).  Therefore,  Gfc  —  1  is  approximately  Gaussian  distributed 
with  zero-mean  and  variance 


Var{Gk  -  1)  = 


( N,N-k ) 

X/  ® n+k,n 


^-I^l)2n=max(1>1.fc) 


(99) 


where  mutual  independence  of  Nn( 0,  l)’s  is  assumed.  Let 


£b  = 


bk 

JVar{Gk-l)' 


(100) 
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we  obtain 


Pr(\Gk  -  1|  >  bk  I  -Ho)  =  2(1  -  *(eb)).  (101) 


If  all  cr^n s  are  the  same,  observing  the  bound  b  in  Theorem  4  completes  the  proof. 
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(c) 

Figure  1:  A  dynamical  nonlinear  neuron  model,  where  neuron  A  is  considered  as  the 
post-synaptic  neuron. 

(a)  Neuron  A  is  influenced  by  presynaptic  neurons  B\,B2,  •  •  • ,  Bn. 

(b)  A  synaptic  connection  between  neurons  A  and  B;  the  influences  of  other  neurons 
on  neuron  A  are  summarized  by  Ut. 

(c)  An  equivalent  probabilistic  version  of  the  neuron  model.  The  impact  of  the 
random  input  Ut  is  now  moved  to  the  spike  generator  where  the  threshold  becomes 
random. 
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Figure  2:  Simulations  for  pair-wise  excitatory  and  inhibitory  correlations. 

(a)  Excitatory  coupling  h(t,s )  =  0.8e-2O(<-s\  t  >  s.  Shown  is  the  two-dimensional 
normalized  scatter  plot  generated  by  the  spike  trains  of  the  two  neurons;  below  it 
is  the  histogram  Gk  that  results  from  collapsing  the  scatter  plot  along  the  principal 

diagonal.  It  corresponds  to  the  function  Nv{k )  =  eb(kh  The  upper  and  lower  bound 
lines  represent  the  95  %  confidence  measure. 

(b)  Inhibitory  coupling,  similar  to  (a)  for  h(t,$)  =  -S.Oe-20^-*),  t  >  s. 
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-64  k=Q  64 


Figure  3:  Interaction  among  three  neurons.  The  network  structure  is  displayed  on 
the  top  graph:  neuron  B  inhibits  neuron  A  and  excites  neuron  C,  and  neuron  C 
excites  neuron  A.  —  —  1.8e_20t,  h^dt)  =  3.6e  20<,  and  hcB(t)  —  2.0e  20t. 

The  top  curve  gives  the  theoretical  connectivity  from  formula  (30)  with  7 *(t,  s ).  The 
middle  one  is  the  correlation  curve  corresponding  to  formula  (52)  generated  from 
spike  trains  A  and  B  only.  The  correlation  is  so  distorted  that  actual  inhibition 
becomes  a  false  excitation  (which  is  actually  due  to  a  strong  excitatory  input  from 
neuron  C).  The  bottom  curve  shows  the  tripartite  correlation  according  to  formula 
(53),  which  displays  the  correct  inhibitory  sign  for  the  connectivity. 
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Figure  4:  Comparison  of  the  preferred  normalization  with  the  difference  normaliza¬ 
tion  (shuffle  method). 

(a)  The  absence  of  a  direct  connection  case  ( h  -  0):  neurons  A  and  B  have  a  common 
input  source  —  a  neuron  driven  by  a  stimulus.  The  connection  strength  from  the 
common  input  is  w  =  1.  The  top  curve  gives  the  collapsed  version  of  the  joint 
PST  histogram  without  any  normalization.  The  correlation  peak  is  purely  due  to 
stimulus  effects.  The  middle  curve  represents  the  difference  normalized  correlation. 
The  bottom  curve  shows  the  preferred  normalized  correlation  curve.  Both  methods 
perform  well  in  indicating  the  absence  of  connection  between  A  and  B. 

(b)  The  presence  of  a  direct  connection  case  (h  f  0):  neurons  A  and  B  have  a 
common  input  source  as  in  (a),  and  in  addition,  a  direct  synaptic  connectivity  from 
B  to  A ,  /mb(0  —  0 .4e~20i.  The  top  curve  gives  the  theoretical  correlation  predicted 
from  Np(k )  =  eh(k\  The  middle  curve  shows  the  difference  normalized  correlation. 
Although  the  connectivity  is  weak  (only  0.4),  the  large  sharp  peak  in  the  correlation 
leads  to  a  false  impression  of  high  excitatory  connectivity,  which  is  in  fact  due  to 
stimulus  effects.  The  bottom  curve  shows  the  preferred  normalized  correlation,  which 
is  very  close  to  the  theoretical  function  ^e_20(- 


